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ABSTRACT 

We investigated the effect of inverse Compton scattering in mildly relativistic 
static and moving plasmas with low optical depth using Monte Carlo simulations, 
and calculated the Sunyaev-Zel'dovich effect in the cosmic background radiation. 
Our semi-analytic method is based on a separation of photon diffusion in 
frequency and real space. We use Monte Carlo simulation to derive the intensity 
and frequency of the scattered photons for a monochromatic incoming radiation. 
The outgoing spectrum is determined by integrating over the spectrum of the 
incoming radiation using the intensity to determine the correct weight. This 
method makes it possible to study the emerging radiation as a function of 
frequency and direction. As a first application we have studied the effects of 
finite optical depth and gas infall on the Sunyaev-Zel'dovich effect (not possible 
with the extended Kompaneets equation) and discuss the parameter range 
in which the Boltzmann equation and its expansions can be used. For high 
temperature clusters (/csT e ^ 15 keV) relativistic corrections based on a fifth 
order expansion of the extended Kompaneets equation seriously underestimate 
the Sunyaev-Zel'dovich effect at high frequencies. The contribution from plasma 
infall is less important for reasonable velocities. We give a convenient analytical 
expression for the dependence of the cross-over frequency on temperature, 
optical depth, and gas infall speed. Optical depth effects are often more 
important than relativistic corrections, and should be taken into account for 
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high-precision work, but are smaller than the typical kinematic effect from 
cluster radial velocities. 

Subject headings: (cosmology:) cosmic microwave background — galaxies: clusters: 
general — methods: numerical — plasmas — scattering 



1. Introduction 

Inverse Compton scattering of the cosmic microwave background radiation (CMBR) by 
hot electrons in the atmospheres of clusters of galaxies, the Sunyaev-Zel'dovich (SZ) effect 
(Sunyaev and Zel'dovich 1980), has become a powerful tool in astrophysics. It is one of the 
most important secondary effects which cause fluctuations in the CMBR. We will refer to 
the effect arising from static gas as the static SZ (SSZ) effect, and that arising from gas 
with bulk motion as the kinematic SZ (KSZ) effect. Fluctuations in the CMBR caused by 
the SZ effects in an ensemble of clusters of galaxies should dominate on angular scales less 
than few arc minutes. The nature of these fluctuations depends on the evolution of clusters, 
and so it is a test of structure formation theories (eg. |Aghanim et al. 1998| ; |Molnar and 
Birkinshaw 1998]). 



Observations of the SSZ effect, begun in the 1970s, have now become routine in the 90s 
with dedicated instruments using the latest receiver technology (for reviews see [Rephacli 



1995b| ; |Birkinshaw 1999]) . SZ effect and x-ray measurements probe the physical conditions 



in the intracluster gas in clusters of galaxies, and allow us to deduce the distance to the 
cluster without additional assumptions. This provides a useful independent method for 
determining the Hubble constant. 

Observations of the KSZ effect are much more difficult since it is typically an order 
of magnitude smaller than the SSZ effect and has the same spectrum as the primordial 
fluctuations in the CMBR. The KSZ effect provides a method of measuring radial peculiar 
velocities of clusters, and even without the tangential velocity component, which might be 
determined using the Rees-Sciama (RS) effect ([Rees and Sciama 1968 ; Birkinshaw and Cull 



1983| |Curvits and Mitrofanov 1986| ; |Aghanim et al. 1998| |Molnar and Birkinshaw 1998| ) 



it should provide important information on large scale velocity fields, which are closely 
related to the large scale density distributions and thus to the average total mass density in 
the Universe. Useful limits on the size of the KSZ effect for two clusters have recently been 
reported by Holzapfel et al. (1997a). 

Most discussions of the SSZ effect have been based on non-relativistic calculations 
of its amplitude, made via a Fokker-Planck type expansion of the Boltzmann equation 
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( [Kompanccts 1957 ). The advantage of this approach is that in interesting cases (electron 



temperature, T e , much greater than the temperature of the incoming radiation) it provides 
a convenient analytical solution for the spectrum of the emerging radiation. However, it 
has been recognized recently that relativistic effects become important for clusters with 
kT e ^ 10 keV. Rephaeli (1995a) provided a relativistic solution for the SSZ effect as a series 
expansion in the optical depth (<^C 1 for clusters). There is no exact analytical solution. 
The numerical integrals involved are tractable in the single-scattering approximation, 
which is usually adequate in clusters. Fargion, Konoplich and Salis (1996) developed exact 
expressions for relativistic inverse Compton scattering of a laser beam with monochromatic 
isotropic radiation, finding good agreement with the approximations of Jones (1968). Their 
expression for the frequency redistribution function (FRDF) for scattering of mono-energetic 
electrons with monochromatic photons agrees with Rephaeli's result. Rephaeli's method 
has been used in a series of papers to evaluate the SSZ effect for hot clusters. Relativistic 
corrections to the SSZ effect, and to the standard thermal bremsstrahlung formulae, were 
applied to determine the Hubble constant in hot clusters by Rephaeli and Yankovich (1997), 
however, their corrections of the thermal bremsstrahlung equation were further corrected 
by Hughes and Birkinshaw (1998). Holzapfel et al. (1997b) used the relativistic results in 
their determination of the Hubble constant from observations of cluster Abell 2163. 

The most general treatment of Compton scattering in static and moving media has 
been derived by Psaltis and Lamb (1997) as a series expansion. As Challinor and Lasenby 
(1998b) noted, however, more terms in the expansion should be taken into account for 
accurate treatment of clusters of galaxies. Recently the Kompaneets equation has been 



extended to contain relativistic corrections to the SSZ and KSZ effects QStebbins 1997 



|Challinor and Lasenby 1998a| , b; Itoh, Kohyama and Nozawa 1998; Nozawa, Itoh and 



Kohyama 1998; [Sazonov and Sunyaev 1998b|) . Starting from the Boltzmann equation, an 
expansion in the small parameters of the dimensionless temperature, = kT e /(m e c 2 ), 
fractional energy change in a scattering, ihv' — hu)/kBT e , and dimensionless radial 
velocity for the KSZ effect, v ra d/c, leads to a Fokker-Planck type equation (the extended 
Kompaneets equation). Corrections up to the fifth order in have been derived (Jltoh 



|et al. 1998 ). These calculations demonstrate the importance of the relativistic effects 



(in accordance with the results of Rephaeli 1995a). Note however, that, as Challinor 
and Lasenby (1998a) emphasized, the extended Kompaneets equation is a result of an 
asymptotic series expansion, therefore it is important to estimate the validity of the 
expansion using other methods. Nozawa et al. compared the convergence of their expansion 
to a direct numerical evaluation of the Boltzmann collision integral, and concluded that 
in the Rayleigh- Jeans region the relativistic corrections give accurate results in the entire 
range of cluster temperatures. Significant deviations are found at higher frequencies for 
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high temperature clusters. 

The SSZ and KSZ effects must be separated in order to extract information on peculiar 
velocities. Fortunately the two effects have different frequency dependence. The maximum 
of the KSZ effect (in thermodynamic temperature units) occurs at about the "cross-over" 
frequency where the SSZ effect changes sign from being a decrement to an increment. In 
a non-relativistic treatment the cross-over frequency is a constant, 218 GHz, independent 
of electron temperature, optical depth, and all other parameters. Rephaeli (1995a) showed 
that in the relativistic case the cross-over frequency depends on the temperature, and his 
results were used by Holzapfel et al. (1997a) in determining peculiar velocities of two 
clusters. Sazonov and Sunyaev (1998b) and Nozawa et al. (1998) give approximations for 
the cross-over frequency as a function of dimensionless temperature and radial peculiar 
velocity. They also conclude that relativistic corrections to the cross-over frequency are 
important, and should be taken into account in future experiments. 

Other methods have been used to investigate inverse Compton scattering, such as 
numerical integration of the collision integral ( |Corman 1970Q , multiple scattering methods 
( |Wright 1979| ), and Monte Carlo simulations. Simulations of inverse Compton scattering 



in relativistic and non-relativistic plasma have been carried out for embedded sources 



(Pozdnyakov, Sobol, and Sunyaev 1983; |Haardt and Maraschi 1993| ; |Hua and Titarchuk 



|1995|) . Gull and Garret (1998) used Monte Carlo methods to evaluate the Boltzmann 
collisional integral. Sazonov and Sunyaev (1998a) used Monte Carlo simulations to derive 
the SZ thermal and kinematic effects. 

In this paper we study the effect of optical depth and non-uniform bulk motion on the 
SZ effect using a Monte Carlo method to calculate the frequency redistribution function. 
The inverse Compton scattering of CMBR photons is treated in the Thomson limit for 
static and infalling plasmas (SSZ and KSZ effects) with spherical symmetry, uniform density 
distribution, and low optical depth over a wide range of gas temperatures and observed 
frequency. We apply our results to clusters of galaxies assuming a static and radially 
infalling (or collapsing) gas component. 



2. The Method 
2.1. Formalism 

The emerging intensity of a beam of radiation in the line of sight after passage through 
a scattering atmosphere can be expressed as a convolution of the FRDF and the incoming 
intensity: 
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I(x) = j F{s)B{v Q )ds, (1) 

where x = hv jk^TcB is the dimensionless frequency, B(vq) is the incoming intensity 
(hereafter assumed to be Planckian), and F(s) is the FRDF, which specifies the probability 
of scattering from Vq to v as a function of the logarithm of the dimensionless frequency, 
u/i/q, s = \n(i//uo), where h, u, ks and Tqb are the Planck constant, the frequency, the 
Boltzmann constant, and the temperature of the CMBR, Tqb — 2.728 ± 0.002 K ( [Fixsen et 
al. 1996|) . The FRDF can be expressed as 

F(s) = e' T 5(s)+w sc P(s), (2) 

where the first term containing the Dirac delta function, S(s), describes the attenuated 
incoming radiation by a factor depending on the line of sight optical depth, r, the 
"out-scattered" radiation, and the second term describes the contribution from scattering 
into the beam, which depends on the FRDF of the scattered radiation, P(s), and a 
weight, w sc , which determines what fraction of the radiation scatters into the beam. This 
decomposition is possible because in our approximation the fractional frequency change is 
independent of the frequency (see equation [16| later). The change of the intensity in the line 
of sight may be expressed as 

AI(x) = (e~ r - l)B(x) + J w sc P(s) B(u )ds. (3) 

For isotropic and homogeneous scattering conditions, so that the scattering parameters do 
not depend on where the scattering happens, the weight 

w sc = (l-e- T ), (4) 

which means that the out-scattered radiation is balanced by the same amount of in-scattered 
radiation, and so there would be no net intensity change (AI(x) = 0) if there were no 
frequency change (P(s) is the Dirac delta function). Our task is to calculate P(s) and 
w sc . However, the assumption in equation ([|) breaks down where the radiation field is 
not isotropic within the cloud, as will be the case in our static and collapsing models, or 
where there is relativistic bulk motion, which introduces anisotropy in the scattering via the 
relativistic beaming effect. In the cases which we discuss in the present paper, equation ([|) 
is an excellent approximation as we have been able to verify using the results of our Monte 
Carlo simulations (see section |2.2| ). Significant departures from equation (f|) will occur 
where the scattering optical depth becomes large, or where the gas velocities approach the 
speed of light: the appropriate treatment in these cases is discussed in a forthcoming paper. 
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Our approximations are adequate for clusters of galaxies, thus we are going to assume the 
validity of equation ([|) in the rest of this paper. 

We derive P(s) using a Monte Carlo method. At low optical depth (r ^ 1), the problem 
is suitable for Monte Carlo simulation because we do not have to follow photons through 
many scatterings in the medium (the average number of scatterings being approximately r) . 
Although in this limit most of the photons do not scatter, and hence provide no information 
on P(s), this is not a problem since we can use the method of forced first scattering (see 
below). 



2.2. Monte Carlo Method 

We give a short description of the method here, for a more detailed description see 
Molnar (1998). 

We assume an isotropic incoming low temperature radiation field (the CMBR). We 
use forced first scatterings to study the inverse Compton process. The photons Compton 
scatter from an electron population with a relativistic Maxwellian distribution of momenta 
in the rest frame of bulk motion. We compute scattering probabilities in the rest frame 
of the electron in the Thomson limit ( Chandrasckhar 1950 ). This involves coordinate 



transformations from the observer's frame to the rest frame of the bulk motion and to the 
rest frame of the electron. We assume time translation invariance and spherical symmetry. 
Time translational invariance is not exact for our model with infall, so that we make a 
snap-shot approximation. The error arising from this approximation is less than the light 
crossing time over the infall time (~ v r /c = j3 r ), which is only a few per cent of the infall 
term for our models. 

In most cases we use the inverse method to generate the desired probability distribution. 
We use a rejection method when the inverse method leads to non-invertible functions or is 
too slow: for a general description of generating probability distributions cf. Pozdnyakov et 
al. (1983); Press et al. (1992). The former reference describes an alternative Monte Carlo 
method to treat inverse Compton scattering. 

In the description that follows we use the word "photon" in the singular to refer to one 
Monte Carlo "photon", one experiment in our simulation. We use a weight, Wi n , to express 
the number of photons this one experiment represents (the weight does not have to be an 
integer). We carried out the simulation in five steps. 



Step 1. The position and direction of incoming photons: 
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We assume that the photons arrive uniformly on a unit sphere (the radius of the gas is 
scaled to unity). In a coordinate system which is placed at the point of impact, the direction 
cosine of the incoming photons from the normal, //j n , can be sampled as 

m n = \f[RN], (5) 

where we use [RN] to indicate a uniformly distributed random number drawn each time 
when it occurs. The azimuthal angle is assumed to be uniformly distributed between and 
2tt. 

Step 2. distance to the forced first scattering: 

We use forced first scattering, which means that we take the probability of scattering equal 
to one on the photons' original line of flight through the gas. Using the inversion method, 
the optical depth at which the incoming photon scatters becomes 

n = - ln(l - [RN] ■ (1 - e- 7 "™-)), (6) 

where r^ a3 , is the maximum optical depth in the line of sight (the superscript refers to the 
number of times the photon has already scattered; zero in this case). Since we are using 
forced first collisions, we have to account for the fraction of photons which are unscattered 
on their path through the cloud. The scattered weight may be obtained from equation @: 
w sc = (1 — e~ T " iax )w in , and stays the same during subsequent (unforced) scatterings. The 
weight of the photons passing through the cloud without scattering is w out = W)„e" T ™ 1 . 

Step 3. scattering: 

At the calculated position of the (forced first) scattering, we use the direction of the 
incoming photon and obtain the direction of propagation and frequency of the scattered 
photon. In the case that the gas is moving, we make a Lorentz transformation into the rest 
frame of the moving plasma. We sample the scattered electron's dimensionless velocity, 
(3 e = v e /c, from a relativistic Maxwellian distribution 

P(P e ) d(3 e = N rel & 7 5 e"^ e d/3 e , (7) 

where the normalization is 

N rel = (eAT 2 (i/e))~\ (8) 

K 2 is the second order modified Bessel function of the second kind, and the dimensionless 
electron temperature is 
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e 



k B T e 



(9) 



We used the rejection method to sample (3 e . 

The distribution of the direction of electron momenta is simplest in a frame in which 
the photon momentum unit vector points into one of the coordinate axes, z for example. In 
this coordinate system the probability distribution of /i e , the cosine of the angle between 
the unit vector of the direction of photon propagation (z axis) and electron velocity, is 



where the sign of the square root was determined so that in the limit of small electron 
velocities we recover the result for an isotropic distribution (/i e = 2[RN] — 1). The angular 
distribution of the scattered electrons in the plane perpendicular to the momentum vector 
of the photon is isotropic (at azimuthal angle uniformly distributed between and 2ir). 
In the rest frame of the electron, the cosine of the polar angle of the incoming photon is 
derived from a Lorentz transformation as 

/i e + j3 e f-tn\ 

" = ww (12) 

where the negative sign in front of fi e is appropriate for an incoming photon. In the 
electron's rest frame the scattering probability of a photon coming in with direction cosine 
ji and leaving with direction cosine /i' is given by Chandrasekhar (1950) 




(10) 



The inversion method leads to sampling /i e as 




(11) 




(13) 



ji' can be sampled using a uniform probability distribution by inversion of 



[RN\ = £' /(//, H'W = 1 ((^ 2 - \)n> Z + (3 - /i V - |) 



(14) 



which leads to a cubic equation for 



(X - + (3 - »y - I (2 [RN] + I) = 0. 



(15) 
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This cubic equation has a single real solution with absolute value of // less or equal to 
one. We now transfer the direction and frequency of the scattered radiation back to the 
observer's frame. The dimensionless outgoing frequency of the photon normalized to the 
incoming frequency (expressed with the s parameter) becomes 



Step 4- l 00 P over scatterings: 

Having the point of scattering, the scattered frequency and the direction of the scattered 
photon, we now sample the optical depth to the next scattering. We do not use forced 
scattering, so the optical depth follows from the usual (inverse) method as 



where the lower index on r refers to the nth scattering (n > 1), and we used [RN] = 1 — [RN], 
which is correct for uniform probability distributions. If r is less than the maximum optical 
depth in the direction of the photon momentum after the previous scattering, r™~^, the 
photon is taken to have scattered within the cloud, and we calculate the new scattering 
direction and frequency of the scattered photon as in step 3. If r n > t™"*, the photon 
escaped (scattered n — 1 times), and we register the impact parameter (the distance between 
the line of sight and the center of the spherically symmetric scattering medium), the weight 
of the photon, the number of scatterings and s, the dimensionless frequency of the escaped 
photon. 

Step 5. The frequency redistribution function: 

At the end of the simulation we sum the weights in every impact parameter bin to 
check our assumption of homogeneous scattering (equation |j) and we determine how the 
average number of scatterings depend on the impact parameter. It turns out that, in our 
case of low optical depth and very mild infall velocities, the average number of scatterings 
has no noticeable dependence on the impact parameter, and that equation is satisfied 
at each impact parameter within the accuracy of our Monte Carlo simulations, which is less 
than 0.1 %. The scattered FRDF depends only on the number of scatterings, thus we can 
sum all the photons to determine an average scattered FRDF, which can be used as an 
excellent approximation to the FRDF corresponding to an arbitrary line of sight. 




(16) 



\n[RN], 



(17) 



Therefore the discrete probability distribution of the scattered FRDF can be derived 
by binning the frequencies of all the out-coming photons as 
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P(s k ) = N s Af k . (18) 

A/fe is the number of photons in the k th bin, for which s k — As/2 < s < s k + As/2, where 
s k is the center of the k th bin, As is the width of the bin, and N s is the normalization, 
N s = 1/ (AsNmc)- Nmc is the number of Monte Carlo photons. This P(s k ) is our sampled 
approximation to P{s), and we then fit an exponential of polynomials to P(s k ) to get a 
convenient expression for the FRDF. Our fit gives an approximation accurate to better 
than half a percent, except in the (small) extended tails, where P(s) is under-represented. 
However, these regions lie 3 orders of magnitude below the peak, and the error arising from 
the fit is negligible for our applications. 



2.3. Testing the code 

2.3.1. Single- scattering approximation 

We tested our code by comparing our single scattering Monte Carlo results for 
the FRDF (Pf IC ) to those derived from Rephaeli (1995a). Rephaeli's single-scattering 
approximation can be written as 

Pi(s) = ~ool\j \]l^K\h + h + fs) d(3 e , (19) 
62JS re i Jf3 

where the lower limit (3q is the minimum f3 e needed to get the frequency shift s, 



M _ i 

A) = nn — -■ (20) 



The functions fi, f 2 and f 3 are 



fi = e 3s [3 3 e (fil - /i 2 - n\ + m) 

h = (ft ~ 3(4e s + 1) + 3e*(:r 2 + Xl ) + ^ ~ 3(1 + ^ )e s (x 2 - Xl ) 

h = 2e s (e s + l)(3-/5 e 2 )ln^, (21) 

Xi 

where x\ — 1 — (3 e /ii, x 2 — 1 — /3 e A*2, A*i and jj, 2 are defined by 



r-i s <o 



- 11 - 



H 2 = < fl> ~ , (22) 

I 1 s > 

and we used equation flltf ) to eliminate This result (equation |19|) agrees with Fargion 
et al. (1996). These expressions can be integrated numerically, except when s = 0. In that 
case Po = 0, and direct numerical integration is not possible because of the diverging /3~ 4 
term. For small j3 e (j3 e < b = 0.1 for example) we can expand the logarithm, and use this 
expansion as a good approximation. For s = 0, equation (|i~9|) becomes 

p ^=v=sh{J° , ° m+ l lf ° d? <)' (23) 

where the integrand is 



fo = 7 e^ /e /?- 4 (2 7 2 Pi + 6p e + (3 - Pi) m(^|) ) . 
A Maclaurin expansion of the logarithm to order ft* gives an adequate approximation 



(24) 



£f M » J\e-^ ( 7 > + I - | + (i - + (i - <„, (25) 

with no remaining divergent terms for the first integral in equation (p3|). 

We derive Pf IC (s) from our Monte Carlo simulation by using the results of only 
the first (forced) scatterings. Figures la and lb show our Monte Carlo results, , 
superimposed on P\ from equations ( |19"1) and (|2l|) for dimensionless temperatures G = 0.03 
and 0.3. The agreement is excellent, confirming that our Monte Carlo code is successfully 
reproducing P\{s). 



2.3.2. Testing the Numerical Integral for the Intensity Change 

We derive the intensity change from the FRDF using a numerical integral (equation |3|). 
The Kompaneets approximation leads to the following FRDF: 

where the Compton y parameter is 



y = J n e a T Od£, (27) 
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where n e is the electron number density as a function of length in the line of sight measured 
by t and G is the dimensionless temperature (for a discussion see for example |Birkinshaw| 
|1999| ; |Molnar 1998| ). In order to check our numerical method, we used the Kompaneets 
FRDF (equation |2BJ) in the numerical integral in equation (|^), and compared the 
resulting intensity change to that of obtained by the analytic solution for the Kompaneets 
approximation 

AI K = y^^( X ^±±-A, (28) 
y (e x - l) 2 V e x - 1 / ( J 

where io = 2(/cb7c-b) 3 / '(he) 2 . We concluded that our numerical method is accurate better 
than 0.1%. 



2.3.3. Relativistic Corrections to Kompaneets equation 

We also compare our results to those from the extended Kompaneets equation up to 
the 5th order in O ( [Itoh et al. 1998[ ). Itoh et al. expressed the intensity change as 



An tnx 3 , , 

A/= 29 

n e x — 1 

and provided expressions for ^ (note, that their y parameter is actually r, the optical 
depth). On Figure |2| we plot AJ/r from the Kompaneets approximation, from Itoh et 
al.'s expansion, and for our single scattering Monte Carlo results. From the figure we 
conclude that our single scattering Monte Carlo result agrees with that of Itoh et al. at low 
temperatures, < 0.03 (T e < 15 keV). Deviations from the Itoh et al.'s result are already 
appearing at = 0.03, and become more pronounced at higher temperatures and high 
frequencies, as we would expect. 

We conclude that our simulation method passes these two tests, and can now be used 
to calculate the effects of multiple scattering and bulk velocity on the SZ effect. 



3. Results 

We performed a number of simulations for uniform density spherical models which are 
either static or have radial infalls with constant gas velocity at all radii. These models 
cover a range of r , the optical depth of zero impact parameter, electron temperature, T e , 
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and gas infall speed, (3 r . Simulations with monochromatic incoming radiation were used to 
determine the scattered FRDF. 

We verified via our simulations that for our low optical depth static models and for our 
models with low optical depth and very mild infall velocities, a parameter space adequate 
for clusters of galaxies, the dependence of the scattered FRDFs on the impact parameter 
is negligible, and that equation (f|) provides a very good approximation to the weight 
of the scattered radiation. Although we determined an averaged scattered FRDF from 
all scattered photons regardless of their impact parameter, in our case, the determined 
scattered FRDF can be used at any impact parameter, since the dependence of the average 
number of scatterings on the impact parameter is negligible (see section |2.2| , Step 5). 

In Figure |3] we show the FRDFs of scattered photons emerging from our spherical static 
models with maximum optical depth r = 0.05 and seven different temperatures. We used 
all photons to derive the scattered FRDFs. At higher temperatures more photons scatter 
into higher energies, thus the FRDFs are broader, and have lower peaks (since they are 
normalized to unity). We show the effect of finite optical depth in Figure f|: higher optical 
depth leads to more scatterings, and therefore more photons scattered to higher energies. 
Even for an optical depth as large as tq — 0.1 the change in the function is relatively small. 
In Figure [| we show the effect on the FRDF of gas infall. Larger infall velocities cause more 
up-scattering of the photons, and hence more spreading of the FRDF, but the most obvious 
change is that the sharp peak at s = is smoothed out by the motion of the plasma. At 
lower temperatures bulk motion causes larger departures from the static FRDF since the 
infall speed is larger relative to the electron thermal velocity. 

We used these results for the scattered FRDF to calculate intensity change using 
equations (|3|) and (^). We evaluated the emergent intensity change at zero impact parameter 
(i.e. through the center of the gas sphere), where r = tq. In Figure ^|we show the intensity 
change AI for a static plasma for two optical depths and five temperatures. Non-zero 
optical depth causes only slight changes in the emerging radiation. Figure [7| shows the 
intensity change for a plasma with infall for two infall velocities and three plausible cluster 
temperatures. Only small changes in the spectrum are apparent, even with such large 
velocities. 

One measure of the spectral deformation that has been used to quantify relative 
correction, and which is of use in determining the frequency at which to search for the KSZ 
effect, is the cross-over frequency. Challinor and Lasenby (1998a) and Birkinshaw (1998) 
suggested a linear expression for how the frequency changes with temperature, as 
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X l in = 3.830(1 + 1.13 6), (30) 



while [Itoh et aL 1998| suggested a quadratic approximation 



X q = 3.830(1 + 1.1674 9 - 0.8533 9 2 ) (31) 

for most cluster temperatures. In Figure || we compare our Monte Carlo results with these 
and other expressions that include relativistic corrections. Our Monte Carlo results for 
single scattering are close to those obtained by numerical integration of the collision integral 
( [Nozawa et al. 1998| ). Relativistic corrections of third and fifth order ( |Challinor and 



[Lasenby 1998a| ; [Nozawa et al. 1998| ), or the linear approximation ( |Challinor and Lasenby 



1998a; Birkinshaw 1999j ) are of varying accuracy in describing the curve: the linear and 



third order expressions give the best results, but extending the series to the fifth order is 
much poorer. This is a consequence of the asymptotic nature of the series, as emphasized 
by Challinor and Lasenby (1998a). 

Figures |9| and [10] show the cross-over frequency as a function of dimensionless 
temperature 9 for finite optical depth and infall velocity. Including a finite optical 
depth causes only a small change in the X (9) curve, and this change does not depend 
much on temperature. Based on our Monte Carlo simulations, we suggest the following 
approximation for the cross-over frequency for single scatterings in static spherical plasma 
for dimensionless temperature 9 £ 0.3: 

X S (9) = 3.827 (1 + 1.2038 9 - 1.2567 9 2 + 0.9098 9 3 ), (32) 
which fits better than 0.01 % in this range with a shift from the optical depth dependence 



AX r (r ) = r (0.35 - 0.04416 9~ a5 ), (33) 

which fits better than 0.005 % for 0.05 < 9 < 0.3, and about 0.1 % for lower temperatures, 
and t < 0.5. 

Figure (^) shows results for our spherical models with gas infall. As for non-zero 
optical depth, additional energy transfers occur because of motion of the gas. As we would 
expect, at low temperatures the contribution to the electron velocity from bulk motion is 
comparable to that from thermal motion, and an enhanced frequency shift results, while 
at high temperatures this contribution becomes negligible. Based on our Monte Carlo 
simulations, we suggest the following approximation for the cross-over frequency shift in 
plasma with infall for 0.01 < 9 < 0.3: 
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AX ^(G,ft) = 0.224 ft 2 G" ' 7 . (34) 
This formula fits the cross over frequency uq better than about half a percent. 

At small optical depth and bulk velocity we can assume that the shifts simply add, so 
that the final expression for the cross-over frequency becomes 

X (G, r , ft) = X*(G) + AXJ(to) + AXo^G, ft). (35) 

In Figure [ll| we show the cross-over frequency as a function of temperature and optical 
depth in the parameter range (G <^ 0.04, r ^ 0.05) important for clusters of galaxies. From 
this figure we may come to the conclusion that, for clusters of galaxies, the optical depth 
effect on the cross-over frequency is more important than non-linear terms in the expansion 
in G. We provide more accurate fitting formulae for this range. From our models we find 



X d (Q, r 0) ft) = 3.829 (1 + 1.1624 G - 0.68948 G 2 ) + 0.14 r + AX ^(Q, ft), (36) 

which fits better than 0.03 % for static models in this range of parameters. The optical 
depth independent first term is similar to the approximation provided by Itoh et al. (1998), 
which was obtained by numerically integrating the collision integral. 

The relativistic corrections to the kinematic SZ effect (cluster radial bulk velocity, v ra d) 
have also been found to be important QNozawa et al. 1998| , |Sazonov and Sunyaev 1998b|) 
The shift was found to be 



Ax kin 



ad I 



300 



Vrad 



a i 



+ 



«2 



e e min 1 (G - 

1-3 



G ) 2 



+ a 3 + a 4 Q + a 5 G z 



(37) 



where v ra d is the radial velocity, G min = 1.654 x 10'*, a\ = 3.857 x 1Q~ 6 , a 2 = —4.631 x 10~ 
a 3 = 1.370 x 10~ 2 , a 4 = 1.014 x 10~ 2 , and a 5 = 0.01 QNozawa et al. 1998; ), which agrees 
well with Sazonov and Sunyaev's (1998b) result (expressed by a different fitting function). 
In Figure [11], long dashed lines represent the results of this kinematic effect of cluster radial 
velocity of ± 100 kms" 1 (Xq 1 (Q,0) + AXQ m (Q,v ra d)) Comparing our results for the shift 
in the cross-over frequency to results from relativistic kinematic effect, we conclude that a 
shift caused by an optical depth of 0.01 is equivalent to a shift caused by a radial velocity 
of about v ra d = — 10 kms -1 . 



We estimate the amplitudes of these effects on the hot cluster, Abell 2163, which 
was discussed by Holzapfel et al. (1997a). The maximum optical depth of the cluster is 
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r = 0.01, the temperature of the intracluster gas is close to 6 = 0.03. Using our results for 
the static effect with the given temperature and maximum optical depth we get about 100 
MHz shift to higher frequencies relative to the linear expression of Challinor and Lasenby 
(1998a). An infall velocity of f3 r = 0.01 causes about an additional 15 MHz shift to higher 
frequencies relative to our result for the static model. These shifts are small relative to the 
20 GHz band width of the instrument of Holzapfel et al. (1997a) and the error in H and 
cluster radial peculiar velocity from ignoring their presence would be about 10 % (if the SZ 
effect is measured at x ^ 5) and 15 kms -1 . By comparison, the component of primordial 
anisotropy in this scale corresponds to adding a velocity noise about ±200 kms -1 . 



4. Conclusions 

We investigated the effect of finite optical depth and bulk motion on inverse Compton 
scatterings in spherically symmetric uniform density mildly relativistic plasma. We assumed 
isotropic incoming radiation (CMBR), a relativistic Maxwellian distribution for the electron 
momenta, and scatterings in the Thomson limit. We demonstrated the usefulness of our 
Monte Carlo method for solving the radiative transfer problem, and calculated the static 
and kinematic SZ effects with different optical depth and gas infall velocities. 

The solution of the extended Kompaneets equation (with corrections up to the fifth 
order) is equivalent to a single-scattering approximation, and significant deviations from it 
occur for hot clusters and at high frequencies. These deviations may be as large as 5 % 
of the intensity change and neglecting them could cause about a 10 % error in the Hubble 
constant. A finite optical depth causes further small changes in the SZ effect: these changes 
may exceed the relativistic correction terms. For typical cluster temperatures, an accurate 
expression for the cross-over frequency as a function of temperature, optical depth, and 
bulk motion is 



As it can seen from Figure |TT|, the cross-over frequency is sensitive to the cluster radial 
velocity, and less sensitive to the finite optical depth. Measurements of the cross-over 
frequency can, in principle, be used to determine the radial velocity of the cluster (e.g., 
as in [Holzapfel et al. 1997b| ), with small extra corrections for optical depth and possible 
gas motion inside the cluster. However, the relatively strong variation of X Q with 0, 
compared to r or /3 r , suggests that the largest uncertainty will arise from the assumption of 
cluster isothermality, even if effects of confusion from primordial (and secondary) CMBR 
fluctuations can be excluded. 



Finally we note that our method can be extended to any geometry, density distribution 
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and complicated bulk motion as desired, and may be used to study the SZ effect in high 
temperature plasmas with or without bulk motion. 

SMM is grateful to Bristol University for a full scholarship, where most of this work was 
done. This work was finished while SMM held a National Research Council-NASA/GSFC 
Research Associateship. We thank to our referee, Dr. Challinor, for suggestions which 
helped us to clarify the presentation of the results. 
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Fig. 1. — The frequency redistribution function, P(s), for a static, single scattering case at 
dimensionless temperatures 6 = 0.03 (a), or 6 = 0.3 (b). The solid line shows the result 
from Rephaeli 1995b) 's semi-analytic method, the boxes with vertical error bars (too small 
to be visible) show results from our Monte Carlo simulations. 
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Fig. 2. — The intensity change AI/r (in units of to = 2(ksTcB) 3 /(he) 2 ) as a function of 
dimensionless frequency x = hv j {\zbTcb) for dimensionless temperatures = 0.01, 0.02, and 
0.03 in static spherically symmetric models. The solid, dashed and dotted lines are results 
from Monte Carlo method (single scattering), from relativistic corrections to the Kompaneets 
approximation (Itoh et al. 1997) and the Kompaneets approximation (Kompaneets 1957) 
respectively. Note that at low temperatures, 6 = 0.01 and 0.02, the single scattering Monte 
Carlo method and the Kompaneets approximation with relativistic corrections give very 
similar results (the solid and dashed lines overlap). 
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Fig. 3. — Scattered frequency redistribution functions, P(s)-s, of scattered photons emerging 
from our spherical static models with maximum optical depth r = 0.05 and seven different 
temperatures. We used all emerging photons to determine P(s) regardless of their impact 
parameter, (see text for details). The higher the temperature, the lower the peak and wider 
the function due to larger energy transfers from the hot electrons to the photons. 
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Fig. 4. — The inner part of the scattered frequency redistribution function, P(s), for static 
clusters with dimensionless temperatures 6 = 0.01, 0.02, and 0.03. The solid and dashed 
lines show the extremes of single scattering and scattering with maximum optical depth 
ro = 0.1. More scatterings lead to more energetic photons, and hence more scattering from 
the line center to the high energy tail. 
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Fig. 5. — The single scattering frequency redistribution function, P(s), for infalling plasma 
with dimensionless temperatures O = 0.01, 0.02 and 0.03. The solid, dashed and dashed dot 
lines show results for models with infall velocities f3 r = (static), 0.03, and 0.05. 
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Fig. 6. — The intensity change AI/t (in units of t = 2(k B T C B) 3 / (he) 2 )) as a function 
of dimensionless frequency x = hv /{k B T CB ) for five dimensionless temperatures in static 
spherically symmetric models, for single scattering (solid lines) and r = 0.1 (dashed lines). 
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Fig. 7. — The intensity change AI/r (in units of i$ = 2(/cbTcb) 3 /(^c) 2 )) as a function of 
dimensionless frequency x = hvj {kB^CB) f° r three dimensionless temperatures in spherically 
symmetric models with infall. Monte Carlo results are plotted using single scattering 
approximation for static gas (solid line) and gas with bulk motion (f3 r = 0.05, dashed line). 
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Fig. 8. — The cross-over frequency as a function of dimensionless temperature, G. The solid 
and long dashed lines show our Monte Carlo results for single scattering, and the results of 
a numerical evaluation of the collision integral (Nozawa et al 1998). The dashed dot line 
represents a linear approximation (Challinor and Lasenby 1998a; Birkinshaw 1998). The 
short dashed and dashed dot dot dot lines use relativistic corrections of third and fifth order 
in 0. 
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Fig. 9. — The cross-over frequency as a function of dimensionless temperature, 6 for static, 
spherically symmetric plasma. The solid and dashed lines are our results for single scattering 
and To = 0.1, respectively. The dashed dot line shows a linear approximation (Challinor 
and Lasenby 1998a; Birkinshaw 1998). Finite optical depths cause a shift due to multiple 
scattering which is almost independent of the temperature. 
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Fig. 10. — The cross-over frequency as a function of dimensionless temperature, 6, for our 
spherically symmetric model with infall. The solid and dashed lines are for infall velocities 
P r = (reference static model) and f3 r = 0.05. We used single scattering results to show 
only bulk motion effects. The dashed dot line shows a linear approximation (Challinor 
and Lasenby 1998a; Birkinshaw 1998). The effect of infall becomes inconsequential at high 
temperature since the infall speed becomes negligible relative to the thermal speed of the 
electron. 
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Fig. 11. — The cross-over frequency as a function of dimensionless temperature, 6, for 
several optical depths in a static spherical plasma. Our Monte Carlo results are from single 
scattering (solid line), r = 0.02 (short dashed line), 0.05 (dash, dot, dot, dot line), and 
0.1 (dotted line). The dash-dot line, which is hardly distinguishable from the solid line, is 
the linear approximation of Challinor and Lasenby (1998a), and Birkinshaw (1998). As a 
comparison, we plot the (large) effect of a radial cluster velocity of ± 100 km s -1 (kinematic 
effect) with long dashed lines (Nozawa et al. 1998; Sazonov and Sunyaev 1998b). 



